% -------------------------------------------------------------------------
function main_dro_
clc;
close all;

% mu = 3.04036e-6; % parameter of system (Sun-Earth)
mu = 1.21506683e-2; % parameter of system (Earth-Moon)
RelTol = 1e-13; % tolerances of integration
AbsTol = 1e-13; % tolerances of integration
MaxNum = 25; % maximum number of iterations
Maxtf = 10; % maximum time of flight

auxdata.mu = mu;
auxdata.RelTol = RelTol;
auxdata.AbsTol = AbsTol;
auxdata.MaxNum = MaxNum;
auxdata.Maxtf = Maxtf;

options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'Events',@event_);

% dro初值
x0 = [-7e4/384400+(1-mu) ; 0 ; 0 ; 0 ; 0.5 ; 0]; %月心


% 打靶法求解dro
[x0_new, period] = dro_(x0 , auxdata); % 微分修正后的初值
[~ , xx0_new] = ode45(@fun_crtbp_ , [0 , 40] , x0_new , options , auxdata); %利用修正后的初值求解crtbp


% 画dro轨道
figure('position' , [100 , 100 , 400 , 280]) % 绘制区域的位置和大小
plot3(xx0_new(: , 1) , xx0_new(: , 2) , xx0_new(: , 3) , 'k' , 'LineWidth' , 2); hold on;
xlabel('x/LU'); ylabel('y/LU'); zlabel('z/LU');
axis equal;
axis([-0.5 2 , -1 1 -0.1 0.1])

% 积分终止条件，即积分到x-z平面终止
function [value,isterminal,direction] = event_(t , xx , auxdata)
value = xx(2); % x-z plane
isterminal = 1;
direction = 1; % direction of event
end

end
